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1. Introduction 

Solution of the six-degree -of -freedom flight equations for aircraft 
and missiles continues to represent one of the most important applica- 
tion areas for analog, hybrid, and digital computer systems. Important 
computer requirements such as accuracy and speed are dependent very 
much on the choice of axis system for the translation equations of motion. 
In this connection it is well known that the flight -jiath axis system 
makes much lower accuracy and speed demands on the computer than 
does the body-axis system [ 1, 2] . Despite this a number of current 
computer mechanizations continue to use body axes for solving the trans- 
lational equations of motion. Because of this, unnecessary demands of 
accuracy or frequency response are placed on the computer and many 
mechanizations which could be all analog or all digital have shifted to 
hybrid implementation. Even if the mechanization is hybrid from the 
outset, there is considerable advantage to be gained by using an efficient 
axis system. The purpose of this paper is to point out again the ad- 
vantages of flight path axes and to summarize the overall equation 
requirements for solving the six-degree of freedom flight equations. 


2 . Body -Axis Translational Equations 


For comparison purposes we present first the body -axis transla- 
tional equations. The body axes x^^, y^^, and are defined as a right- 
hand set fixed to the vehicle with the Xj^ axis along the longitudinal axis 
and the axis directed downward for normal level flight. The com- 
ponents of the total vehicle velocity vector along the Xj^, y^^ and 
axes, respectively, are Uj^, Vj^, and (see Fig. 2.1), The components 
of the body -axis angular velocity vector n (and hence the vehicle angular 
velocity vector) along the y^^, and axes are and i. e. , 

roll rate, pitch rate, and yaw rate, respectively. Here we assume 
Vp and n represent vehicle translational and rotational velocity vectors 
as viewed from an inertial (non -accelerating) frame of reference. If 
we denote the external forces along a set of coordinate axes by X, Y, 
and Z, respectively, thenEuler’s translational equations of motion, ob- 
tained by summing forces along the coordinate axes, are the following: 

m(U-VR+WQ)-X (2.1) 

m(V-WP+UR)-Y (2.2) 

m(W-UQ + VP) = Z (2.3) 

where m is the mass of the vehicle. The inefficiency of these equations 
in body axes is immediately apparent when one considers the approximate 
size of the various terms. Let the vehicle be, say, a Mach 2 aircraft 
with = 2000 ft/ sec. A reasonable upper limit on pitch -rate Q 


2 




might be 2 radians/ sec. Thus, the term in Eq. (2.3) might 

2 

be as large as 4000 ft/ sec or 125 g’s ! On the other hand Z^/m, 

b ’ 

the normal acceleration due to the external force (primarily gravity 
and aerodynamic lift) may have an upper limit of several g's. Hence 
artificial accelerations which are perhaps 20 to 50 times greater than 
the actual accelerations are introduced because of the high rotation 
rates which the body -axes experience. This means much less favorable 
computer scaling and hence much poorer solution accuracy for a given 
computer precision. Fiirthermore, the high-speed dynamics of the 
rotational equations are coupled into the translational equations, thus 
placing severe dynamic response requirements on the computer. The 
use of flight -path axes greatly alleviates these problems. As we shall 
see in the next section, the flight path axes allow a more efficient 
calculation of the aerodynamic angle of attack a and the aerodynamic 
angle of sideslip ^ than body axes allow. Using body axes, and assuming 
that the ambient air mass is not moving relative to the inertial frame 
used to define V , then the following formulas can be used to obtain a , 

Yr 

and velocity magnitude V from the body-axis velocity components 

Jr' 

Wb 

tan a = — (2.4) 
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(2.5) 


1 o b 
sin^ = ;^ 

P 

/ 2 2 

V = 5 = (u^ + V, + W, ) 

p cos a cos /3 \ b b b / 

3. Flight -Path Axis Translational Equations 

Next consider the flight -path axes x , y , z shown in Fig. 3. 1. 

w w w 

These differ from the body axes Xj^, yj^, z^ by the angle of attack a and 
the angle of sideslip as shown in the figure. To rotate from body 
axes to flight -path axes one first pitches the body axes about y^ through 
-a. This defines an intermediate axis system x , y , z called 

S S o 

stability axes. One then yaws about z , through /3, which defines the 

s 

fl^ht-path axes x , y , and z . Note that V is the x,„ component of 
vehicle velocity; the y and z components are zero by definition. 

AV TV 

Let us define the x , y , and z components of flight -path angular 
w w w 

velocity relative to inertial space by P^, Q^, and and the com- 
ponents of external force along x , y , and z by X , Y , and Z , 
respectively. Then, since U = V and V = W = 0, the translational 

AV P tV Aw 

equations (2. 1), (2. 2) and (2. 3) referred to the flight -path axes become 


mV = X 

(3.1) 

p w 

mV R = Y 

(3.2) 

p w w 

mV O = Z 
p ^ w 

(3.3) 
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Solution of these three equations results in total velocity V^, flight - 

path axis yaw rate R^, and flight -path axis pitch rate Q^. 

Next consider the formulas for a and /3, Reference to Fig. 3. 1 

shows that a is directed along y with a component a cos ^ along y . 

s 

W 

Thus a cos /3 is equal to the difference between body-axis and flight -path 
axis angular rates along y . Therefore, from Eq. (3. 3) we can write 


o cos ^ = Q cos /3 - P ® sin /3 + — ^ 
m b ^ mV 


(3.4) 


where P, is the body -axis (not stability axis) angular rate along x 
” E 

and is given by 


P, = P. cos a + R. sin a 
b b b 


(3.5) 


Similarly, reference to Fig. 3. 1 shows that ^ is directed along 
and is equal to the difference between flight -path axis and body axis 
angular rates along z . Thus from Eq. (3.2) 


w s 

mV ^b 
P 


(3.6) 


where R^^ is the body -axis angular rate along z and is given by 


R^ ~ ® ^ ^b ^ 


(3.7) 


(^ote that 
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Figure 3. 1 Flight Path Axes x^, y^, and the 
Relation to Body Axes x^^, y^^, 
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(3.8) 



since the y. body axis and y stabilityaxes are coincident. ) 

Equations (3. 1), (3. 4) and (3. 5) can be integrated to yield total 
velocity V^, angle of attack a, and angle of sideslip jS. They do not 
present the scaling difficulty of Eqs. (2. 1), (2. 2), and (2.3) in the body 
axes. The body-axis velocity components U, , V^, and can be ob- 

DO D 

tained from V^, o', and /3 by the formulas 


U, = V cos a cos B 

bp 

(3.9) 

V, = V sin /3 
b P ^ 

(3.10) 

W. = V sin a cos B 
bp ^ 

(3.11) 


Similarly, the flight -path -axis forces X , Y , and Z can be derived 
from body -axis force components Yj^, and by the formulas 

^s " ^ ^ ^b = Xg cos ^ + Yg sin (3.12) 

^s=^b ’ Y^ = -X^sin^ + Y_ cos^ (3.13) 

^s ° + ^b'°® ® ’ ^w'^s 

where X , Y , and Z are the intermediate stability -axis force com- 
ponents. Frequently the aerodynamic force components are computed 
along stability axes, in which case the power plant and gravity forces, 
computed in body axes, would be resolved into stability axes where the 
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aerodynamic forces are added; then the total forces would be resolved 

into flight -path axes to allow use of Eqs. (3. 1), (3. 4), and (3. 5). It 

should be noted that we have assumed throughout this section that the 

translational and rotational velocity vectors are velocities relative to an 

inertial reference frame. If the atmosphere through which the vehicle is 

flying can be considered to be fixed with respect to this inertial frame, 

then the velocity magnitude is the vehicle velocity relative to the 

atmosphere, and a and /3 represent the aerodynamic angle of attack and 

sideslip, respectively. Using the approximation that the earth is flat, 

and with constant surface winds, it is possible to define the inertial 

reference frame as a frame attached to the atmosphere. Then all the 

formulas, as presented, are correct and a, ^ and can be used for 

computation of aerodynamic forces and moments. Unfortunately, for a 

rotating spherical earth axes fixed in the atmosphere are not inertial. 

If we consider such a frame to be inertial, we will make acceleration 

A2 , A 

errors in Eqs. (3. 1), (3.2) and (3.3) of the order of V /r^ where V is 

the vehicle velocity relative to an inertial frame with origin at the center 

of the earth and r is the radius of the earth. For illustration, consider 
0 

a vehicle flying eastward in still air with a velocity V relative to the 
atmosphere. Then V “= + r^cOj^cos L, where is the earth spin 

rate and L is the latitude. If the vehicle is flying westward, 

A ^1 

V = V - r a; , cos L. For other headings V lies between these values, 
a o N 


9 



For example, if V = 3000 ft/ sec and L = 0 degrees, V ranges between 

A 

approximately 1600 and 4400 ft/ sec. The corresponding acceleration 

A2 

error in Eqs. (3.1), (3. 2) and (3.3), given by V /r^, ranges between 

2 

approximately 0. 1 and 1 ft/ sec . For many flight vehicles this is a 
negligible error. On the other hand for a supersonic transport cruising 
eastward at 3000 ft/ sec this lowers the required steady -state lift by 
about 3 per cent, which could lower the drag significantly and hence 
make a noticeable difference in maximum range. 

There appears to the authors to be no simple way to take these 
accelerations into account and still use a flight -path axis system re- 
ferenced to the atmosphere. One could add an approximate correction 

2 

acceleration in the vertical direction given by V / r to the translatory 

P ® 

forces in Eqs. (3. 1), (3.2), and (3.3). In fact, one could further sim- 
plify the computation by adding the term only to Eq. (3, 3) based on the 
argument that most of the time a supersonic aircraft will be in near- 
level flight at cruise, and that a moderate acceleration error dviring 
transient conditions can be accepted. 

In this case Eq. (3. 3) becomes 



Note that this equation will exhibit acceleration errors of the order of 
2 , 

V Vr for conditions far from level flight, 
p 0 ^ 
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4. Rotational Equations 


The only reasonable axes to use for the rotational equations of 
motion are the body axes. If we let the external moment components 
along x^, y^j and be Mj^, and Nj^, respectively, then summation 
of moments about the three body axes of a body symmetrical about the 
Xj^ plane leads to the equations : 

■xx - Oyy - W «b ''b - 'xz " Pb S,> = ^ 

lyy - «zz - W «b ^b " 'xz f^b' - = “b 

I R, - a - I ) Pu Ok - I , (Pix - Oh = N, (4.3) 

zz b XX yy b ^ xz b ^ b b 

Here 1,1, and I are the moments of inertia about x, , y, , and 2 l 
XX yy zz b b b 

respectively, and I is the product of inertia of the symmetrical body. 

xz 

Note in Eqs. (4. 1), (4. 2), and (4. 3) that the second term in each equation 

represents a nonlinear inertial coupling term. For flight vehicles such 

as large transport aircraft which do not generate relatively high angular 

rates these terms often can be neglected. For many flight vehicles the 

roll -rate P. has a maximum value which is considerably higher than 
b 

pitch-rate or yaw -rate Rj^. Hence the term in Eq. (4. 1) often 

can be neglected in comparison with the ^.nd Pj^j^ terms in 

Eqs. (4. 2) and (4.3), respectively. 

The third term in each of Eqs. (4. 1), (4. 2), and (4. 3) represents the 

effect of the proauct of inertia I . If the x, y, z body axes have been 

xz 
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chosen to be almost coincident with the principle axes, this term may be 

negligible in all three equations, since I will be very small compared 

with the principle moments of inertia. For relatively low angular rates 

2 2 

the nonlinear terms “ ^b ’ S)^b^ usually can be neg- 

2 

lected and in any event frequently can be neglected as small compared 
with in Eq. (4. 2). 

5. Computation of Euler Angles for a Flat Earth 

Solution of Eqs. (4. 1), (4, 2) and (4. 3) results in computation of the 

body-axis angular velocity components Pj^, and R^^. These must be 

integrated a second time to obtain the orientation of the vehicle body 

axes with respect to the desired references axes, typically Euler axes 

which point North, East, and toward the center of the earth (x , y , 

e -^e 

and in Fig. 5.1). This orientation usually is expressed in terms of 
the conventional aircraft Euler angles, i. e. , heading -angle ijy, pitch 
angle 9, and bank angle 0. These angles usually are computed from 
^b’ ^ ^b following well-known equations: 


^ - (R^ cos <j> + sin 0)/ cos 6 

(5.1) 

0 = ^ cos 0 - R|^ sin 0 

(5.2) 

0 = Pj^ + 0^ sin 6 

(5.3) 


Note that, since R^^ are the body axis components of the 
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vehicle angular velocity relative to an inertial reference system, there 
is a small error introduced by the angular velocity of axes which point 
North, East, and down relative to a spherical earth. This small error 
can be corrected, if necessary, using Eqs. (7. l)-(7. 3) in the next 
section. 

The well known singularity of the Euler angle system at 6 = ±n / 2 
can be avoided, if necessary, by computing direction cosines or 
quaternions [4] instead of Euler angles, or by introducing a fourth 
angle [ 3] . 

The use of quaternions rather than direction cosines should be 
considered if a system free of singularities is needed, since there are 
only four quaternions with a single redundancy as compared with nine 
direction cosines with six redundancies. 


6. Computation of Vehicle Position for a Flat Earth 

Once the orientation of the flight vehicle with respect to the Euler 
axes has been established, e.g. , by means of the Euler angles, then 
it is possible to compute the velocity north, U ; the velocity east, V , 

G 0 


and the velocity downward, W^, or its negative, the rate of climb h. 
Direct integration then yields the vehi^*le position. 

Determination of U^, V^, is complicated by the fact that the 


vehicle velocity vector, V lies along the x wind axis and therefore must 
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be resolved from wind to earth axes. Unfortunately, the complete 

orientation of the wind axes is known only relative to body axes, hence 

it is necessary to perform the resolution of Vp into earth axes by first 

resolving it into body axes, then from body axes to earth axes. The 

resolution of V from wind axes to body axes is accomplished by the 
P 

transformation given in Eqs. (3. 9), (3. 10), and (3. 11). 

The resolution of vehicle velocity from body axes to earth axes can 

be accomplished by using direction cosines. Thus let f fgj ^3 

be the projections of a unit vector along the x body axis onto the x^, 

y^, and earth axes, respectively. Similarly, let m^, m^, and m^ 

be the projections of a unit vector along the y body axis onto the x^, 

y , and z earth axes, respectively. In the same way let n., n^, and 

Hg be the projections of a unit vector the z body axis along the x^, y^, 

and z earth axes, respectively. Then by definition 
e 

velocity north “ Uj^ + m^^ + n^ (6. 1) 

velocity east = ^2 ^b ^ ^2 ^b ^ ^2 ^b 

velocity downward = ~ ^ “ ^3 ^ *^3 ^b ^ ^^3 ^b 

It is easy to show that the direction cosines are related to Euler angles 
by the following formulas: 
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= cos 9 cos J// 


= cos 0 sim// I (6. 4) 

Xg = - sin 9 ' 

nij^ = - cos (j> sin x}y + sin 0 sin 9 cos p 
nig = cos 0 cos J// + sin 0 sin 0 sin i// ( (6. 5) 

mg = sin 0 cos 9 

n^ = sin 0 sin 0/ + cos 0 sin 9 cos 4' 
ng = - sin 0 cos 4^ + cos <f> sin 9 sin 4^ 
ng = cos 0 cos 0 

(Equivalent formulas for direction cosines in terms of quaternions are 
given in Ref. 4. ) An alternative mechanization of Eqs. (6. 1) through 
(6. 6) avoids computation of the direction cosines by instead performing 
successive resolution of the velocity components Uj^, Vj^, and 
through the angles - 0, - 0, and - 4^. Consider the intermediate axis 
system x', y’, z' in Fig. 5. 1. Clearly the velocity components U’, V, 
and W along x', y', and z' are given by 

U' = Ub , V = cos 0 - sin 0 , W' = sin 0 + cos 0 (6.7) 

Defining U”, V", and W" as the velocity components along the inter- 
mediate axes x”, y”, and z”, we have 
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U" = U' cos 0 + W’ sin 0 , V" = V’ , W" = - U’ sin 0 + W’ cos 0 (6.8) 


Finally, from Fig. 5. 1 we see that 

U = U" cos - V” sin , V = U” sin i// + V" cos ,W = -h=W” (6.9) 

Q B ^ 

Successive application of Eqs. (6.7), (6.8), and (6.9) for U, V, and 

W to obtain U , V , and W requires fewer mathematical operations 
e e e 

than using Eqs. (8. 1) through (8. 6). It therefore has computational 
advantages using either an analog or digital mechanization. In com- 
puting ground coordinate.s from U and V in Eqs. (6. 1) and (6,2), or 

G G 

Eq. (6. 9), it is important to note that U and V represent airspeed com- 

ponents north and east, respectively. To convert them to groundspeed 

components the north component of wind, w , must be subtracted from 

U , and the cast component of wind, w , must be subtracted from V , 
e y ^ 

respectively. Thus if s and s represent distance traveled north and 

X y 

east, respectively, then 


s 


X 




( 6 . 10 ) 



( 6 . 11 ) 


Equations (6. 10) and (6. 11) are valid only for steady winds, since the 


implicit assumption was made that axes stationary relative to the atmos- 
phere are inertial. Equations (3. 1), (3.2), and (3. 3) are referred to 
inertial space and correction terms must be added if the reference axes 
are not inertial. 


17 



7. Computation of Vehicle Euler Angles and Position for a Rotating 
Spherical Earth — — __ 

In the previous section we presented the formulas for computing 

vehicle position over a flat earth with steady winds. We can use the 

same position formulas to obtain velocity north, U^, and velocity east, 

V , over a rotating spherical earth with radius r , angular rate co . 
e ON 

However, U and V will represent airspeed components and must be 
corrected to yield groundspeed components s^ and s^ respectively. 
Furthermore, when a spherical earth is considered it may be necessary 
to correct the vehicle angular rates used to compute Euler angles in 
order to take into account the rotating reference frame, as pointed out 
earlier in Section 6. It can be shown that the body axis components 
of the vehicle angular velocity relative to the Euler reference frame, 
*’be> «be- «be are given by the following formulas [ 5] : 



s s 


s 

+ -^ tan L 
r o 


(7.1) 





-^1 
r 1 


+ -U2+^tanLn3 


(7.3) 


Here r is the radial distance of the vehicle from the center of the earth 
and is given by 
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r = r + h 
o 


(7.4) 


where r is the radius of the earth and h is the vehicle altitude. In 
o 

many cases we can substitute r^ for r in Eqs. (7. 1), (7.2), and (7.3) 

and still obtain sufficient accuracy. 

The values of P,. , Q , and R, in these equations are then used to 
be ^e be 

compute the Euler angle rates. Thus by analogy with Eqs. (5. 1), (5. 2), 
and (5.4) 

cos <f) + sin (f>)/ cos 6 (7.5) 

e = cos <t> - Rj^^ sin 0 (7 . 6) 

i = P, + ^ sin 9 (7.7) 

In many cases the computations involved in Eqs. (7. 1), (7.2), and (7.3) 
can be neglected, i. e. , we can assume that ~P^^, 0,^^ = Qj^, 

This is particularly true if the overall six-degree -of -freedom computa- 
tion involves a control system (automatic or human) which attempts to 
maintain 6, and at specified values. In any case the correction 
rates are the order of Vp/ r^. For example, if Vp = 2000 ft/ sec, the 
correction rate is equal to approximately 0. 005 degrees per second. 

On the other hand, if the flight -vehicle problem includes a stable plat- 
form, the rate corrections given by Eqs. (7, 1), (7.2), and (7.3) may be 
important. 
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It should be noted that s and s in Eqs. (7. 1), (7,2), and (7.3) 

X y 

represent vehicle velocity components north and east, respectively, 
over the surface of a non-rotating earth with steady winds. On the other 
hand Eqs. (6. 1) and (6.2), or alternatively, Eqs. (6,9), gives us 
and Vg, i. e. , vehicle velocity components north and east, respectively, 
relative to the inertial reference frame for the translational equations 
of motion. We made the approximation in Section 3 that this reference 
frame is fixed relative to the ambient atmosphere. We can account 
for the linear velocity of the atmospheric reference frame (but not the 
angular velocity) by noting that relative to the surface of a non-rotating 
earth it is moving northward with the northerly component of wind and 
eastward with the sum of the rate due to earth spin and that due to east- 
ward component of wind. Thus we can write the following equations: 


s = U - w 
X e X 

(7.8) 

s = V + 03,r r cos L - w 
y e N y 

(7.9) 

It can also be shown [ 5] that the time rate of change of latitude and 

longitude are given by the following formulas: 


s 


L = — 

r 

(7.10) 

s 


^ T -^KT 

r cos L N 

(7.11) 


20 



It should be noted that motion of a flight -vehicle over a rotating 
earth can be treated exactly [ 5] , but that the exact translational 
equations referred to axes fixed relative to the ambient atmosphere, 
are very complicated. Thus many of the computer mechanization ad- 
vantages for flight-path axes are lost, and the computation of ff , /3, and 

V is much less elegant. We have attempted in this section to describe 

P 

how one can utilize the flight -path axis system and still correct approxi- 
mately for the fact that the vehicle is flying over a rotating earth with 
surface winds. This approach should be adequate for all but the most 
exacting requirements unless the vehicle reaches hypersonic speeds. 

For subsonic vehicles or supersonic vehicles traveling over relatively 
short distances the flat -earth equations in Sections 3 through 6 should 
be adequate. 

8. Computation of Aerodynamic Forces and Moments 

Aerodynamic forces and moments for flight vehicles normally are 

computed in stability axes. Thus let X^, Y^, and be the aerodynamic 

forces along the x , y , and z stability axes shown in Fig. 3. 1. Here 
s s s 

- X corresponds to drag D and - Z corresponds to lift L. 
a. ^ 

In addition to aerodynamic forces the flight vehicle experiences 

propulsion and gravity forces. These are most conveniently defined 

in body axes. Let us denote the propulsion force components along the 

X and z body axes by X and Z , respectively. The gravity force 

P P 
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components along x, y, and z will by definition be mgX , mgm„, and 
mgHg respectively. The sum of the propulsive forces and gravity forces 
along body axes can then be resolved to stability axes, where the aero- 
dynamic forces are added to obtain the total force component X , Y , 

s s 

and Zg along the stability axes. Using (6. 4), (6. 5), and (6. 6) to express 
the direction cosines f g, m^, n^ in terms of Euler angles, we obtain: 

Xg = (Xp - mg sin 9) cos a + (Z^ + mg cos 9 cos <l>) sin a - D (8. 1) 
Yg = mg cos 9 sin ^ + (8. 2) 

Zg = - (Xp - mg sin 9) sin a + (Z^ + mg cos 9 cos <(>) cos a - L (8. 3) 

Finally, reference to Fig. 3. 1 shows that the force components X , Y , 

w w 

and Z^ along the flight -path axes can be computed using the following 
formulas : 


X = X cos /3 + Y sin /3 
w s s 

(8.4) 

Y = - X sin /3 + Y cos 
w s s 

(8.5) 

z = z 

w s 

(8.6) 


These force components are used to mechanize the translational equations 
given in Section 3. If the aerodynamic forces are given in body axes 
rather than stability axes, the modification of Eqs. (8, 1), (8.2), and (8. 3) 
is apparent. 
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The external moments acting on the flight vehicle consist of aero- 
dynamic moments L , M , and N , normaUy given in stability axes, 

and power plant moments L , M , and N , given in body axes. The 

P P P 

total moments L, M, and N in body axes become the following: 


L = L 

cos a - N sin a 

+ L 

(8.7) 

a 

a 

P 


M = M 

+ M 


(8.8) 

a 

P 



N = L 

sin a + N cos ct 


(8.9) 

a 

a 

P 



Again, if the aerodynamic moments are given in body axes, the sim- 
plification of Eqs. (8.7), (8.8), and (8.9) is obvious (set a = 0). 

Figure 8. 1 shows a block diagram of the over-all six -degree -of - 
freedom equations for the case of a flat earth. 
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Figure 8. 1 Block Diagram of Combined Flight ’Path Axis, Body ’Axis System for a Flat Earth, Steady Winds. 
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